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Abstract 

We present a maximum-likelihood analysis for estimating the angular distribution of power in 
an anisotropic stochastic gravitational-wave background using ground-based laser interferometers. 
The standard isotropic and gravitational-wave radiometer searches (optimal for point sources) are 
recovered as special limiting cases. The angular distribution can be decomposed with respect to 
any set of basis functions on the sky, and the single-baseline, cross-correlation analysis is easily 
extended to a network of three or more detectors — that is, to multiple baselines. A spherical 
harmonic decomposition, which provides maximum-likelihood estimates of the multipole moments 
of the gravitational-wave sky, is described in detail. We also discuss: (i) the covariance matrix of 
the estimators and its relationship to the detector response of a network of interferometers, (ii) a 
singular-value decomposition method for regularizing the deconvolution of the detector response 
from the measured sky map, (iii) the expected increase in sensitivity obtained by including multiple 
baselines, and (iv) the numerical results of this method when applied to simulated data consisting 
of both point-like and diffuse sources. Comparisions between this general method and the standard 
isotropic and radiometer searches are given throughout, to make contact with the existing literature 
on stochastic background searches. 
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I. INTRODUCTION 



Data from the laser interferometric gravitational-wave detectors LIGO [1, 2, 3, 4], Virgo 
[5, 6], and GEO [7, 8] are currently being analysed for the presence of gravitational waves 
from a variety of sources. These include signals from inspiraling and coalescing compact 
binaries (for example, neutron stars and/or stellar mass black holes) [9, 10, 11, 12], contin- 
uous gravitational waves from quasi-periodic sources such as pulsars [13, 14, 15, 16], and 
bursts of gravitational radiation associated with gamma-ray bursts [17, 18, 19], core-collapse 
supernovae, or other violent events [20]. In addition, searches are ongoing for the presence 
of a background of stochastic gravitational radiation of either astrophysical or cosmological 
origin [21, 22, 23], whose detection might provide insights about the very early universe [24], 
well before the production of the cosmic microwave background [25]. 

Although no direct detections of gravitational waves have been made to date, the most 
recent data taken are of unprecented sensitivity [13, 17, 26], leading to upper limits on 
gravitational-wave strengths that are competitive with or surpass those from electromagnetic 
or particle physics observations. Of particular note is the upper limit on the strength of a 
gravitational-wave signal from the Crab pulsar [27], which is a factor of 1.6 lower than 
the corresponding limit inferred from electromagnetic pulsar spin-down observations [28]. 
Also, the current direct limit on the strength of an isotropic stochastic gravitational-wave 
background at 100 Hz Q gw < 6.9 x 10~ 6 [26] (at 95% confidence) has surpassed bounds set by 
considerations from Big Bang Nucleosynthesis [21] and from the microwave background [29]. 

In this paper, we describe an analysis method that estimates the angular distribution 
of power in an anisotropic stochastic gravitational-wave background. This method includes 
both the standard isotropic [23, 30, 31] and gravitational-wave radiometer searches [21, 32] 
(optimal for point sources) as special limiting cases. (For our purposes anisotropic is taken to 
mean not necessarily isotropic.) Similar to the radiometer technique, the method presented 
here looks for modulations in the gravitational-wave signal induced by the Earth's rotational 
motion relative to an anisotropic background. The method provides maximum-likelihood 
estimates of the angular distribution of gravitational- wave power V(Q) = ^ Q P a e Q (f2), 
decomposed with respect to some set of basis functions on the sky. By choosing a pixel 
basis e^,(f2) = 5(Cl,Cl'), we recover the results of the radiometer method discussed in [21, 
32]. By choosing the spherical harmonics basis Yi m (Cl) defined with respect to the Earth's 
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rotational axis, we obtain maximum-likelihood estimates of the multipole moments Vi m of 
the gravitational-wave sky. This basis is particularly convenient as the standard isotropic 
analysis corresponds to simply restricting attention to the monopole moment Pocb while the 
point-source radiometer results are well-approximated by choosing a sufficiently large value 
of /max (^max ~ 30), appropriate for the diffraction-limited beam pattern at / ~ 1 kHz. In 
addition, the use of spherical harmonics simplifies the problem of removing the 'smearing' 
effects of the beam pattern from the measured sky map (that is, deconvolution of the dirty 
map), given the smaller number of elements and symmetries of the beam pattern matrix 
with respect to the Im indices. The problem of deconvolving a cross-correlated gravitational- 
wave signal from the interferometers' beam pattern in the spherical harmonic basis has been 
described in [33]. We address this problem in detail in section IV. We further note that 
the spherical harmonic basis is useful for the efficient analysis of cross-correlated data in a 
variety of applications including searches for transient gravitational- wave sources [34]. 

Regardless of basis, the method described here is easily extended to work with a network 
or three or more detectors with uncorrelated detector noise, by simply adding the individual 
baseline beam patterns and dirty maps before deconvolution. A multi-baseline analysis 
improves the overall sensitivity of the search by reducing the variances of the individual 
estimators, and provides a natural way of regularising the deconvolution of the dirty map; 
the beam pattern matrix has fewer null (or nearly null) directions for multiple baselines and 
is thus more stable during inversion. 

The structure of the rest of the paper is the following: In section II, we briefly review 
the statistical properties of an anisotropic background, and show how a generalized overlap 
reduction function arises in a cross-correlation search for such a background. In section III 
we derive the optimal estimators of the angular distribution of the gravitational-wave power, 
starting from the likelihood function for cross-correlated data. We explicitly construct the 
beam pattern matrix, and discuss its relation to the covariance matrix of the estimated 
V a . Section IV describes details of the data analysis implementation and issues related to 
deconvolution and regularisation. It also briefly describes how to extend the analysis to 
a network of three or more detectors, and the expected increase in sensitivity from using 
multiple baselines. In section V we present numerical results of the method applied to 
simulated data. We consider both point-like and diffuse-source injections, and compare the 
extracted and injected sky-maps. Finally, in section VI we summarize our results. We 
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also include three appendices: Appendices A and B contain definitions of the spherical 
harmonics and some useful identities relating different multipole moments, beam pattern 
matrix components, etc.; Appendix C defines a related detection statistic that assumes a 
particular distribution of (normalized) angular distribution functions on the sky. 

II. GRAVITATIONAL- WAVE BACKGROUNDS 

Stochastic gravitational-waves are produced by the superposition of a large number of 
weak, independent, unresolved gravitational- wave sources. The signal can be either cos- 
mological or astrophysical in nature, leading to different expected characteristics: (i) A 
cosmological background, consisting e.g., of remnant gravitational waves left over from the 
very early universe, is expected to be predominantly isotropic, similar to that of the 2.73 K 
temperature distribution of the cosmic microwave background, (ii) An astrophysical back- 
ground, on the other hand, produced by more recent astrophysical events, such as early-phase 
compact binary inspiral or continuous radiation from pulsars (see, e.g., [35]), will most likely 
be anisotropic, following the spatial distribution of the sources. In addition, cosmological 
backgrounds are expected to have relatively smooth, monotonic power spectra (for example, 
falling off close to f~ 3 for standard inflationary models) [36], while astrophysical backgrounds 
are expected to have power peaked at some characteristic frequency. Not surprisingly, to op- 
timally search for these different signals requires different search algorithms, adapted for the 
angular distribution and spectral properties of the source. In this section, we describe how 
the anistropy of a stochastic gravitational-wave background manifests itself in the statistical 
properties of the signal, and in the expected value of the cross-correlation of the output of 
two detectors. The following sections then describe how one can search for such a signature 
in the measured data. 

A. Statistical properties 

In the transverse-traceless gauge, the metric perturbations due to a stochastic 
gravitational-wave background can be written as a superposition of plane waves having 
frequency / and propagating in the direction fi: 




(2.1) 
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where e^ b (Cl) are the gravitational- wave polarization tensors. (Summation over polarization 
indices A is understood.) In standard angular coordinates on the two-sphere 9 e [0,n], 
(f) G [0, 27r), we can write 

Q = sin#cos0x + sin#sin0?/ + cos 6* 5 , (2-2) 
I = cos 9 cos (fix + cos# sin (fiy — sin 9 z , (2.3) 
rh = — sin (fi x + cos (fiy , (2.4) 

so that {l,m,Q} forms a right-handed system of unit vectors. We can then define the two 
(A = +, x) polarization tensors to be 

e ab(&) = l ah ~ m a m b , (2.5) 
e ^(^) = Lrhb + rhJb ■ (2.6) 

Note that there is a rotational degree of freedom in the definition of polarization tensors as 
one is free to rotate I and m by an angle ip in the plane orthogonal to Cl. For a gravitational- 
wave source with a symmetry axis, such as an inspiralling binary, the angle i/j can be in- 
terpreted as the polarization angle of the source. However, as we will assume that the 
stochastic background is unpolarized, there is no loss of generality in taking ip = 0, so that 
the polarization tensors have the form given above. 

The Fourier coefficients h,A{f, ^) are complex functions that satisfy h^i—f, £1) = h A (f, Cl), 
since h a b(t,x) is real. For a stochastic gravitational- wave background these coefficients are 
random fields whose expectation values define the statistical properties of the background. 
Without loss of generality we can assume that the fields have zero mean: 

(M/,n)) = o. (2.7) 

We will also assume that the background is unpolarized, Gaussian, and stationary, but allow 
for an anisotropic distribution. The most general form of the quadratic expectation value 
satisfying these requirements is 

(h* A (f, n)M/', fr)> = \hl n) s(f - f')5 AA ,5(n, ft) , (2.8) 

where V(f, Cl) specifies both the spectral and angular distribution of the background. The 
factor of 1/4 has been included so that for an isotropic background H(f) = V(\f\,fl) is the 
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one-sided strain power, when summed over both polarizations. Given the above definitions, 
it follows that 

<W/) = -% = %if I dnv(f,ti), (2.9) 

where dp gw is the energy density contained in the frequency interval df. Here Hq is Hubble's 
constant, and p c = 3c 2 Hq/8ttG is the critical energy density needed to close the universe. 
(To prove Eq. (2.9), one should write p^ in terms of an expectation value of the product 
of the time derivatives of the metric perturbations h a b(t,x), and then expand the metric 
perturbations in terms of the plane wave components as in Eq. (2.1), using Eq. (2.8) to 
evaluate the expectation value; see, e.g., [37, 38].) Thus, the energy density in a stochastic 
gravitational-wave background has contributions from all parts of the sky as encoded in the 
all-sky integral of V(f, Cl). 

In what follows, we will assume that V(f, Cl) can be factorized into a product of two 
functions 

V(f,{l)=V({l)H(f), (2.10) 

where 8(f) is a dimensionless function of frequency, normalized so that H(f R ) = 1, where 
Jr is a reference frequency, typically taken to equal 100 Hz (a frequency in LIGO's most 
sensitive band). V(ti) specifies the angular distribution of gravitational- wave power, and 
H(f) its spectral shape. This factorization does not amount to a loss of generality if one 
restricts attention to small enough frequency bands. For our analysis, we will assume that 

H(f) = (f/f R y, (2.ii) 

where (3 is a power-law index which we fix (for example, (3 = for constant strain power). 
Using Eqs. (2.9) and (2.10), one can show that this assumption for H is consistent with 

n gw (f) = n R (f/h) 3+ ^ (2.12) 

where Q R is the fractional energy density in gravitational waves evaluated at the reference 
frequency f R . 

The angular distribution function V(Cl) can be expanded in terms of a set of basis func- 
tions on the two-sphere according to 

V{9) = V a e Q ({l) , (2.13) 
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where summation (or integration) over a is understood, and 



v a = / dnv(n)e* a (n) , (2.14) 

l S 2 



Jaf3 



[ dtt e* a {Cl)e p {Cl) . (2.15) 
Js 2 



f S 2 

The choice of basis, in principle, should not affect the physical search results. However, in 
practice, such a choice can bear on computational costs of a search and also on the systematic 
errors affecting observations results, e.g., arising from the truncation order of the spherical- 
harmonic basis. For these reasons, we expect that while searching for gravitational-wave 
point sources, a decomposition with respect to a pixel basis 

V(Q) =V Cl ,5(tt,tt') (2.16) 

is the natural choice. For a diffuse background, e.g., dominated by a dipole or quadrupolar 
distribution, a spherical harmonic decomposition may be the better choice: 

V(Cl) = V lm Y lm (Q), (2.17) 



v lm = / dnvityYfjn) , (2.18) 

Js 2 

where the second equality follows from our normalization convention for the Y[ m (see Ap- 
pendix A). Note that the pixel basis coefficients, Va, have units of strain 2 /Hz whereas 
the coefficients in the spherical harmonics basis, Vi m , have units of strain 2 /Hz/rad. This 
normalization convention also implies 



27T 2 , 

n* = -^pjl ■ (2.19) 

Note that only the monopole moment Poo contributes to Qr (and hence to figw(/)) as ai l 
higher-order multipole moments give zero when integrated over the sky. 



B. Overlap factor 

We will denote the time-series output of two detectors / = 1, 2 by 

si(t) = h I (t)+n I (t), (2.20) 
where rij(i) is the detector noise and hi(t) is its response to a gravitational- wave background: 

/OO /' 
df / dnh A (f,n)F I A (n,t)e t2nnt - Cl ^ t)/c) . (2.21) 
■00 Js 2 
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Here 

Ff(£l,t)=df(t)e*(fl) (2.22) 

is the detector response function, which encodes the directional sensitivity of detector / to a 
plane-polarized gravitational wave propagating in direction Cl and xj specifies the location of 
interferometer /. (The absolute value \Ff(Q,t)\ plotted as function of direction Cl is called 
the detector antenna pattern.) The detector tensor is 

df(t) = \ [xf{t)XKt) - Y?{t)Y?{t)\ , (2-23) 

where Xj(t), Yi(t) are unit vectors pointing along the interferometer arms for detector /. 
The vectors xj(t), Xj(t), and Yj(t) are all time- dependent due to the Earth's rotation. (We 
are using equatorial coordinates, with the spatial origin at the center of the Earth, f-axis 
pointing along the Earth's rotation axis, and x-axis pointing in the direction of the vernal 
equinox.) 

Given a time-series sj{t), we define its short-term Fourier transform §i(f,t) by 

~ Sl (f,t) = f +T/2 dt' e-^ ft ' Sl (t') , (2.24) 

Jt-r/2 

where r is much greater than the light-travel time between any pair of detectors, but is 
small enough that the the detector response function Ff(Cl,t) and detector location xj{t) 
do not vary significantly with time over the interval [t — r/2,t + t/2]. Typical values of r 
are from a few tens of seconds to a few hundred seconds. The cross-correlation between the 
output of the two detectors is then defined in terms of these short Fourier transforms as 

C(f,t) = -sl(f,t)s 2 (f : t). (2.25) 

T 

The factor of 2 is a convention consistent with the definition of one-sided power spectra, so 
that the total cross-power for a particular time t is given by integrating C(f, t) over positive 
frequencies. These cross-spectra are the starting point for the maximum-likelihood analysis 
described in the following section. 

If the noise at the two detectors is uncorrelated — a reasonable assumption for spatially- 
separated detectors — then it follows that the expectation value of the cross-spectra depends 
only on the gravitational- wave signal components 

(C(f,t))=*(hl(f,t)h 2 (f,t)). (2.26) 



Using Eqs. (2.8), (2.10), and the short-term Fourier transform of (2.21), one can then show 
that 

(C{f,t)) = H{f) f dfi 7 (6,/,t)7>(6), (2.27) 
Js 2 



where 



7(6, /, t) = -Ff(Cl, t)Ff(n, f) e M-(a(t)-a(*))/« . (2.28) 



The function 7(0, /, t) is a geometric factor that takes into account the separation and 
relative orientation of the two detectors (see e.g., [39]). For an isotropic background, 
(C{f,t)) ex Htf)>y{f), where 

T (/) = A/" dQF 1 A (fi,t)F 2 A (fi,t)e i27r/fi - (5l(t) - 52( * ))/c (2.29) 

87T ^2 

is the standard overlap reduction function [40, 41]. The factor of 5/87T is a normalization 
constant chosen so that j(f) = 1 for all frequencies for a pair of coincident and coaligned 
interferometers with 90-degree opening angle between the interferometer arms. Note that 
7(/) is time- independent due to the all-sky integration. 

For subsequent analysis, it will be convenient to rewrite the right-hand side of Eq. (2.27) 
in terms of an integral over the components of V(Cl) and f,t) with respect to a set of 
basis functions on the two-sphere 

(C(f,t)) = H(f)<y a {f,t)V a , (2.30) 

where 

la (f,t) = [ dfi7(^/,*)ea(fi), (2.31) 
Js 2 

j{Cl,f,t) = 7 «(/,*)e*(n). (2.32) 
Therefore, in the pixel basis, a <-> Cl', one has 

7 (Cl,f,t)= 1Cil {f,t)6{Cl,fl'), (2.33) 
and in the spherical harmonics basis, a <-> Im, one gets 

J,t) = 7 / m (/,t)^(fi), (2.34) 

7im(f,t) = [ dttj(nj,t)Y lm (tt). (2.35) 
Js 2 
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Note that the above definition of jim(f>t) differs from that (2.18) of Vi m by a complex 
conjugation, but it agrees with the convention used in [38]. The time-dependence of 7z m (/, t) 
is particularly simple: 



In addition, for an isotropic background, the above definitions imply 7(/) = (5/ V4n)joo(f, t) 
for any t. 

III. MAXIMUM-LIKELIHOOD ANALYSIS 

In this section, we derive the maximum-likelihood estimators of the angular distribution 
V a of the stochastic gravitational-wave power. The analysis given below also makes clear 
the relationship between the covariance matrix of these estimators and the beam pattern 
matrix for the cross- correlation measurements. Since maximizing the likelihood is equivalent 
to minimizing the squared deviation of the estimators away from their expected values, the 
estimators so obtained are optimal in the sense of maximizing the expected signal-to-noise 
ratio of the estimators. Hence, the likelihood analysis reproduces the standard results of 
optimal filtering for the isotropic (e.g., [23, 30, 31]) and radiometer [21, 32] analyses without 
explicitly introducing a filter function Q in the construction of a statistic. 

A. Maximum-likelihood estimators 

The maximum likelihood analysis for an anisotropic stochastic background takes as its 
fundamental data vector the cross-spectra 



evaluated at a set of discrete times t and discrete (positive and negative) frequencies /. As 
shown in the previous subsection, the expectation values of the cross-spectra are given by 





c ft = c(f,t) = -r 1 {f,t)s 2 (f,t) 




(C ft ) = H(f) la (f,t)V a . 



(3.2) 



The covariance matrix is given by 



N 



ftj'f - 



(Cf t Cf, t ,) - (Cf t )(Cf, t ,) 
6w6ff>PiU,t)P2{f,t) 



(3.3) 
(3.4) 
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where Pi(f,t), 1—1,2 are the one-sided power spectra of the detector output for time 
segment t, which satisfy 

<W.*)W,0> = ls ttl 5 ffl Pj(f,t) (3.5) 
« (nfotfaif,*)). (3-6) 

We have assumed that there is no cross-correlated noise, and that the cross-correlated and 
auto-correlated gravitational-wave signal power are much less than the detector noise power 
to obtain the approximate relations Eqs. (3.4) and (3.6). 

Treating Pj(f,t) and the gravitational- wave spectral shape 8(f) as known quantities, 
the likelihood function is then 

1 



p({C /t }|{7> a })ocexp 
where 



^(V) 



(3.7) 



X 2 (V) = (°}t ~ (C} t ))Nj t ] f ,,(C rt , - (C rt ,)) . (3.8) 
tft'f 

Using Eqs. (3.2) and (3.4), we have 



x 2 cp) = *) - H(fmf, t)r: 



i 



-(C(f,t)-8(f)^(f,t)Vp). (3.9) 



P 1 (f,t)P 2 (f,t) 

Since maximizing the likelihood with respect to V a is equivalent to minimizing chi-squared, 
one can show that the maximum likelihood estimators for the V a are given by 

Va=(r- 1 )af)Xp, (3.10) 

where 

^EEff') fl| g fil) ^'). p.") 

(We will adress the invertability of Y a p in section IV.) Note that the standard estimator of 

the strength of an isotropic stochastic background [23, 30, 31] 

-l 

Yt 

at I " at 
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has the same form as the above, with ^ f Y t /o^ playing the role of Xp and l/of the role 

of r a/3 . 

For later reference, we note that the minimum value of chi-squared is xLm = X 2 iJ > )i 
which can be written explicitly as 

2 VV \C(f,t)\* 

- v* a x a - x;v p + Kv aP Vp . (3.14) 

Also, in analogy with what is done for cosmic microwave background experiments such as 
WMAP [42], we can construct estimators of the angular "power" spectra 

1 1 

Cl = aTPT E lVlm]2 (3 - 15) 

m=— I 

by simply replacing Vi m with the estimators V a evaluated in the spherical harmonics basis — 
i.e., 

i ' 

m=—l 

Note that the C\ defined above are actually estimators of the squared angular power (since 
Vim already has units of power), unlike the cosmic microwave background data for which 
the Ci really do have units of power. Also, we will see in the next subsection that these 
estimators are biased. Unbiased estimators of the C\ are given in equation (3.24). 



B. Error estimates 

If the spectral shape H(f) that we assumed for the signal model exactly matches that of 
the observed background, it is fairly easy to show that the estimators V a constructed above 
provide unbiased estimates of the angular distribution of gravitational- wave power: 

(V a )=V a . (3.17) 

This follows immediately from the fact that 

(X a ) = r af) Vf}, (3.18) 

which in turn shows that the X a are the components of the so-called 'dirty' map — i.e., 
X a represent the gravitational-wave power on the sky as seen through the beam matrix 

13 



of the two detectors, T a p. Equation (3.10) shows that by inverting r a p, one obtains the 
components of the 'clean' map, P a . The process of going from the dirty map to the clean 
map is an example of deconvolution. 

In addition, one can show in the weak-signal approximation that 

(x a x;)-(x a )(x;)*r aP , (3.19) 
(^)-(^>(^>«(r-%. (3.20) 

Thus, r a/ 3 is the covariance matrix of the dirty map X a , and {T~ l ) a p is the covariance matrix 
of the clean map V a . A matrix like T a p, whose inverse is the covariance matrix of the signal 
parameters, is often called a Fisher information matrix. An alternative definition of T a p, 
illustrating its connection to the likelihood function, is 

d 2 \np({C ft }\{V a )} 



a/3 



(3.21) 



dV*dVf3 

As is evident from the above expression, if one has several independent measurements (so 
that the combined likelihood is just a product of individual likelihoods), the T a p matrices 
simply add. 

Finally, using the above expressions for the expectation value and covariances of the V a , 
one can show that 

(Ci) « Cj + ^E^" 1 )^. ( 3 - 22 ) 

m 

(Cf) - {Q} 2 » — Kr- 1 )^! 2 • (3.23) 

m,m' 

Note, in particular, that the estimators C\ are biased. Unbiased estimators of the C\ are 
given by 

^ = ^-2iTiX)( r ~ 1 )^- ^ 

m 

C. Decomposition in terms of pixel basis or spherical harmonics 

The analysis presented above has been written in terms of the components V a and 7 a (/, t) 
of V(Cl) and 7(f2, /, t) with respect to an arbitrary set of basis functions on the two-sphere. 
For most purposes, we will be interested in the components with respect to only two bases: 
The pixel basis, for which a <-> Cl and and 7q(/, t) are given by (2.16) and (2.33), and 
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the spherical harmonics basis, for which a <-> Im and Vi m and ji m (f,t) are given by (2.18) 
and (2.35). Each basis has its own set of advantages and disadvantages, which we briefly 
describe below. 

In the pixel basis, is an estimate of the true gravitational-wave power Vq coming 
from direction Q. It is a real quantity and should be non-negative. The quantity X$, on 
the other hand, is the power coming from direction Cl as seen by the detector. It includes 
gravitational-wave power from other directions on the sky due to the finite acceptance of 
the beam pattern function, as well as from instrumental noise. The matrix T^, connects 
the two via (3.18), and can be directly interpreted as a point spread function. It specifies 
how a point source at Cl is spread to other points Cl' by the response of a pair of detectors. 

In the spherical harmonics basis, the Vi m are estimates of the true multipole moments 
Vim of the gravitational- wave power on the sky. The matrix ri m ji m i, is no longer directly 
interpretable in terms of a point spread function, but it plays an analogous role as the 
inverse matrix {T~ l )i m ,Vm specifies the correlations between the various multipole moment 
estimates. 

In addition, the Fisher matrix Y a p has two symmetries: parity (see Eq. A7) and rotational 
symmetry around the z-axis. Since spherical harmonics respect these symmetries, this leads 
to some simplifications. Parity is an exact symmetry, because the only difference between 
gravitational-wave signals coming from antipodes is an opposite sign of the time shift between 
detectors. Therefore the detector noise, as expressed by the Fisher matrix, is identical for 
antipodes. This implies that Ti m ii m i = for all odd I — /' (almost half of the matrix 
elements). Z-axis rotational symmetry is broken by daily variations in detector sensitivity, 
but still implies that Ti m ^ m r m for m ^ mf, i.e., Ti m ^ m i is a block-diagonally dominant 
matrix. The pixel basis has no such symmetry. 

Furthermore, in the spherical harmonics basis it is simple to specify a resolution cut-off 
by only allowing / < Z max . This avoids over-sampling and reduces the number of required 
basis vectors. Also, since extending this cut-off to a larger / max does not affect the original 
basis vectors, it is straightforward to run the analysis with a higher resolution, and later do 
the matrix inversion at a lower resolution. 

Finally, the computationally dominant part of the analysis is the calculation of the Fisher 
matrix. Since the Fisher matrix has N 2 elements, with N the number of basis vectors, 
working in the spherical harmonics basis makes the analysis significantly more efficient. 
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And the mentioned symmetries help to reduce the computational load even more. 



IV. IMPLEMENTATION AND ANALYSIS DETAILS 

As shown in Sec. Ill A, the maximum-likelihood estimators of the angular distribution of 
power in an anisotropic gravitational-wave background are given by 

V a =(T- 1 ) afi Xp, (4.1) 

where Xp are the components of the 'dirty' map (3.11), and (r -1 )^ are the components 
of the inverse of the beam pattern matrix T a p (3.12). In this section, we describe: (i) 
some of the implementation details related to the calculation of Xp and T a p, (ii) a method 
for regularizing the inversion of T a p, (iii) how to extend the single baseline analysis to a 
network of detectors, and (iv) a Bayesian model selection scheme for determining Z max for 
the spherical harmonic decomposition. For concreteness we consider a network of detectors 
consisting of the LIGO interferometers HI and LI and the Virgo interferometer, VI. 



A. Calculating Xp and T a p 

The components of the 'dirty' map Xp and the beam pattern matrix T a p are given by 

^EEff') fl(f S f ,») ^"' < 4 - 2 ' 
r- = EE^*u§fe •»</.«)■ («) 

These are the fundamental data products of this analysis from which (r _1 ) Q/ g and V a are 
then calculated (4.1). Although the various quantities entering Xp and T a p have already 
been defined in the previous two sections, we describe here in more detail how they are 
calculated in practice. 

(i) H(f) = (///r) is the assumed spectral shape of the gravitational-wave background 
(2.11). This is an input to the data analysis pipeline which we fix at the start of the analysis. 
The parameters /r and j3 are the reference frequency and spectral index for the (assumed) 
power-law behavior of the gravitational-wave spectrum. For the analyses described later 
in this paper, we choose /r = 100 Hz and j3 — 0, corresponding to constant strain power. 

1(3 



Other values of fn and (3 are, of course, possible. For example, (3 = — 3 corresponds to 
constant fractional energy density f2 gw (/) = const, which follows from Eq. (2.12). 

(ii) C (f, t) are the cross-spectra of the data, calculated as a product of the short-term 
Fourier transforms of the time-series output of the two detectors, cf. Eqs. (2.25), (2.24). 
The time-series data are first downsampled to a few kilohertz (from 16384 Hz to 2048 Hz for 
LIGO; from 10000 Hz to 2000 Hz for Virgo), high-pass filtered above 40 Hz (to reduce con- 
tamination from low-frequency seismic noise), and then windowed (to avoid spectral leakage 
of strong instrumental lines), before being discrete- Fourier-transformed to the frequency do- 
main. As r is typically of order 100 s (r = 60 s for the simulations that we will describe in 
Sec. V), the frequency resolution of si(f,t) and C(f,t) is of order 1/r = 0.01 Hz, which is 
much finer than what is needed for the other frequency- series H(f), Pi(f,t), and 7 Q (/, t), 
which are typically taken to have a frequency resolution Af = 0.25 Hz. Hence, to match Af, 
we average together several frequency bins of C(f,t). This averaging or "coarse graining" 
has also been used for previous stochastic searches, see e.g. [30]. It is a technique used to 
avoid unnecessary frequency resolution, especially in Pj(f,t) and 7«(/, t). 

(iii) Pj(f,t) are the power spectra associated with the individual (/ = 1,2) detector 
outputs (3.5). We use Welch's modified periodogram method to estimate the power spectra, 
averaging together periodograms from 4-sec long, 50% overlapping, Hann-windowed data, 
which are taken from the two time segments immediately preceding and following — but not 
including — the analyis segment. (The 4-sec data stretches corresponds to the Af = 0.25 Hz 
frequency resolution mentioned earlier.) This technique greatly reduces a bias that would 
otherwise result from non-zero covariance between Pj(f,t) and C(f,t). 

For an actual analysis of real data, one needs to consider the possibility of short-term 
variations in the detector noise that are not consistent between the analysis segment and 
the two neighboring segments from which the Pi(f,t) were estimated. For the simulations 
that we will describe in Sec. V, the data were stationary, so no consistency cut needed to 
be applied. 

(iv) j a (f,t) are the components of the overlap factor j(Cl,f,t) defined by Eqs. (2.28), 
(2.33), (2.35). These are geometric factors that encode the relative separation and orien- 
tation of the two detectors, as specified by the detector response functions Ff(£l,t) (2.22), 
(2.23). In the pixel basis, the components 7q;(J, t) can be efficiently calculated by using 
one Fast Fourier Transform, and reading out the resulting cross-correlation time series at 
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the time shift corresponding to each pixel [21, 43]. In the spherical harmonics basis, the 
components 7/ m (/>^) can be efficiently calculated using analytic expressions derived in [38]. 
In particular, the authors in [38] show that for sidereal time t — 0, one can write 7z m (/, 0) as 
a simple linear combination of spherical bessel functions j n {x)/x n (for I even) or j n (x) / 'x n_1 
(for I odd), where x depends on the relative separation of the detectors x = 2ttJ\xi — X2.\/c. 
The coefficients of the linear combinations are complex numbers that depend on the relative 
orientation of the detectors. Explicit expressions for a few of the 7z m (x) = 7j m (/, 0) for the 
LIGO Hanford-Livingston pair are given below: 



l00 (x) = -0.0766j (x) - 2.1528- 



x 



+ 2.4407^^ , (4.4) 

x 2 



ll0 (x) = -0.0608iji(a;) - 2.6982z j2 



x 



+ 7.72172^^, (4.5) 

x 1 



7ll ( x ) = -(0.0519 + 0.0652z) ji(x) 

. jo (x) 

- 1.8622 + 1.0516z 

x 

+ (4.0106 - 2.4936i) 3 ^$~ ■ ( 4 - 6 ) 

x l 

(Note that the numerical coefficients above do not agree with those in [38] , due to an overall 
normalisation by 47r/5 and phase factor e l7n ^, where (j) — —38.52° is the angle between the 
separation vector between the LIGO Hanford and Livingston detectors and the Greenwich 
meridian.) For arbitrary sidereal times t, one uses Eq. (2.36), which follows from the e im ^ 
dependence of the spherical harmonics Yi m (Ct) = Yi m (9,<f>). Here (9,<f>) are related to the 
equatorial coordinates (ra,dec) via 9 = tt/2 — it (dec/180 ) and <fi = ir (ra/12 hr). 



B. Deconvolution and regularization 

Equation (4.1) is a formal description for estimating the angular structure V a of a 
gravitational-wave background. We refer to this as deconvolution since it tries to remove the 
smoothing introduced by the point spread function. Deconvolution requires inverting the 
Fisher matrix T a p. However, in practice, the Fisher matrix T a p is somewhat ill-conditioned. 
There are two reasons for this. 
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First, the detector pair is diffraction limited. Thus, as we choose a basis with higher 
spatial resolution, the condition number of the Fisher matrix T a p gets worse, resulting in 
a reduced signal-to-noise ratio for the deconvolved map. We can address this by picking 
a basis with a reasonable resolution cut-off, which makes the spherical harmonics basis set 
with / < Z max a natural candidate. 

Second, there are certain power distributions V a to which the detector pair is essentially 
blind. For those distributions positive and negative contributions from different sky locations 
to the total cross-correlation essentially cancel. Mathematically they are described by 

X Q = V aP V p « 0, (4.7) 

i.e., they are the eigenfunctions of the Fisher matrix Y a p with the smallest eigenvalues. 
These eigenfunctions tend to have z-axis rotational symmetry because the detector pair 
is rotating with the Earth. However this symmetry can be broken by daily variations in 
detector sensitivity. To address this second issue, we have chosen to use a singular value 
decomposition (SVD) regularization scheme, which we describe in some detail below. 
Since T a p is Hermitian, its SVD has the form 

T = USU* , (4.8) 

where U is a unitary matrix and S = diag(sj) is a diagonal matrix with non-negative entries 
Si (the eigenvalues of T). Without loss of generality we can further assume that the diagonal 
elements of S are sorted from the largest to smallest values. Then the problematic modes 
according to Eq. (4.7) correspond to the last entries of the diagonal of S. Figure 1 shows the 
relative size of the eigenvalues of a typical T a p matrix in the spherical harmonics basis (with 
l max = 20), taken from the no-injection simulation of Section V. We can now set a threshold 
s m i n on the size of the eigenvalues, setting all eigenvalues Sj < s m i n to infinity (their inverse 
to zero), or alternatively to Sj = s m j n . Using this modified matrix S' we can then define the 
regularized V as 

r' = US'U* (4.9) 

and its inverse as 

r'- 1 = us'- x u* . (4.io) 

The threshold s min is chosen by weighting the quality of the deconvolution (larger point 
spread function for higher values of s min ) against the addition of noise due to poorly measured 
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modes (lower values of s m j n ). While one can make this trade-off argument more quantitative, 
the choice will be somewhat influenced by the spatial shapes one is looking for. For the 
purpose of this paper, we simply chose to keep 2/3 of all eigenmodes, and set all the small 
eigenvalues equal to s m i n . This is somewhat arbitrary, but a reasonable choice to get rid of 
the extremely small eigenvalues of a typical Fisher matrix (Figure 1). As can be seen in 
Sec. V, this choice allows for a reasonable recovery of simulated injections. 

Using this regularization scheme has two side effects that need to be mentioned. First, 
Eq. (3.10) is replaced by 

P a = {T'- l ) aP Xp. (4.11) 
Thus the expectation value of V' a is 

{P a ) = (V- 1 ) aP r^P J ^V a . (4.12) 

This constitutes a bias in the estimator, which is expected since we chose to ignore the 
modes of V a that are poorly measured. Under the assumption that we know the shape of 
the source this bias can be calculated. Assuming the signal consists of point sources, Figure 
2 shows the size of that bias as a function of sky position. Second, the covariance matrix of 
V'a (in the weak-signal approximation) is now given by 

(P'>;> - (V'aHV'p) = (r'- 1 ) a7 r 7 ,(r'- 1 )^ . (4.13) 

Finally, we note that adding additional detector pairs with different baselines can, to a 
certain degree, act as a natural regulator, simply because one detector network might be 
more sensitive to a particular mode than another as illustrated in Figure 1. This is described 
in more detail in the following subsection. 

C. Multiple baselines 

As shown explicitly in [32] for the case of the directed radiometer method, the above 
analysis can easily be extended to a network of three or more detectors with uncorrelated 
detector noise. One simply adds the dirty maps X T a J and Fisher matrices T 1 ^ for each 
distinct detector pair I J: 

i j>i i j>i 
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FIG. 1: Eigenvalues of typical Fisher matrices r Q/ g for different baselines and the multibaseline 
detector network. For this analysis / m ax = 20, corresponding (Z max 

+ l) 2 = 441 total modes. 

For each individual baseline some of the SVD eigenmodes are (almost) null [see Sec. IV B]. The 
multibaseline network, however, has fewer null modes, illustrating the fact that a network of 
detectors acts as a natural regularizer - independent baselines tend to complement each other. 
The plot was produced using the simulated data described in Sec. V. 

| 3 

do.? 
Id is 

FIG. 2: Magnitude of the bias due to the SVD regularization scheme, for the case of point sources. 
A value of implies that the expectation value of the corresponding pixel is equal to the point 
source signal strength, while 1 implies that a point source at that location would not be seen. 

where the subscript M signifies a network of baselines. This follows from extending the 
likelihood formulation in Sec. Ill to include sums over baselines as well as frequency and 
time. The maximum-likelihood estimators V a then retain the same form as for the single 
baseline case, namely 

V a = [(T^UXff. (4.15) 
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FIG. 3: Standard deviation for spherical harmonics components, without SVD regularization. It 
illustrates how the multiple baselines (solid line) reduces the estimation error by natural regular- 
ization. 

This follows immediately from Eq. 3.21. 

Different baselines in the network partly complement each other and help fill gaps in 
sensitivity present in individual baselines pairs. This has an important consequence. The 
sensitivity gaps correspond to degeneracies in the Fisher information matrix, which make it 
hard to estimate the true stochastic background. By filling these gaps the network acts as 
a natural regularizer, as illustrated in Fig. 1. 

Our next step is to quantify how the natural network regularization and the SVD reg- 
ularization reduce the error in estimating the spherical harmonic components. The full 
covariance matrix of the estimated multipoles Vi m is the best measure of estimation error, 
but it is inconvenient to compare covariance matrices. Rather, we use the standard deviation 
of each multipole as our figure of merit for estimation error: 



for regularized estimators (no summation over Im in either of these two formulas). We plot 
the standard deviations for each multipole for both unregularized and regularized estimators 
in Fig. 3 and Fig. 4, respectively. Both of these figures indicate that multiple baselines 




(4.16) 



for unregularized estimators, and 




(4.17) 
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FIG. 4: Standard deviation for spherical harmonics components, with the SVD regularization 
described in Sec. IV B. While the combination of network and SVD regularization leads to the 
minimum estimation error, the SVD regularization alone can significantly reduce the estimation 
error for a single baseline. 

vastly reduce the estimation error. In addition, Fig. 4 illustrates that though the estimation 
error is minimized when both the network and SVD regularization are present, the SVD 
regularization alone can significantly reduce the estimation error (to just ~ 25% more than 
the regularized network error) for the LIGO only baseline. Thus, even if one interferometer 
is not unusable for some period, the regularized spherical harmonic moment estimators for 
the remaining baseline can still provide reasonable results. 

D. Model selection for spherical harmonic decomposition 

In addition to choosing the cutoff for the SVD regularization of the Fisher matrix (as 
described in Sec. IV B), one needs to specify the value of / max , the maximum value of 
the spherical harmonic index I used in the spherical harmonic decomposition. Choosing 
/ max fixes the total number of multipole moments Vi m , and hence defines the signal model. 
Larger values of Z max mean finer angular resolution of the sky maps and more parameters 
available to fit the data. But since the estimators Vi m are correlated with one another, 
increasing the number of parameters simultaneously increases the uncertainty associated 
with each parameter. Thus, there is a tradeoff between accurately modeling the data (more 
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parameters) and minimizing uncertainties (fewer parameters). In this subsection we outline 
how Bayesian model selection can be used to fix Z max . (This discussion is meant to motivate 
future study as we do not implement a model selection scheme in this work.) 

Bayesian model selection (see, e.g., [44]) is a framework in which the data themselves 
determine which signal model is most appropriate. The basic idea is to compare the various 
models (e.g., Mi and M 2 ) by computing the ratio of the probability of the models given the 
data D. By Bayes' theorem, we have 



where p(Mi) and p(M 2 ) are the a priori probabilities of the models, and p(D\Mi) and 
p(D\M 2 ) are the likelihood functions for the data given the two models. The ratio of the 
likelihoods p(D\Mi)/p(D\M 2 ) is known as the "Bayes factor" (see, e.g., [44]). If there 
is no a priori reason to prefer one model over the other (as is often the case), then 
p(Mi)/p(M 2 ) = 1, implying that the posterior odds is just the ratio of the likelihood func- 
tions, p(D\Mi)/p(D\M 2 ). Since a given model often involves a set of parameters a, calcu- 
lating the likelihood of the data for a given model requires marginalizing over the possible 
values of these parameters — i.e., 



where p(a\M) is the prior probability distribution of the parameters for that model. 

In situations where the data is informative — i.e., when the likelihood function p(D\a, M) 
is peaked relative to the prior p(a\M) — we have the approximate relation 



where a is the value of a that maximizes the likelihood, 5a is the range of parameter values 
over which the likelihood is peaked, and Aa is the full range of parameter values. The factor 
5a/ Aa penalizes a model that uses more parameter space volume than needed to fit the 
data. This factor can be understood in terms of Occam's razor, which says that everything 
else being equal, simpler models that can adequately fit the data are preferred. 

For example, if we ignore the subtleties described in Sec. IV B related to the inversion of 
the Fisher matrix Ti m ^ m i, an / max = 30 map (961 parameters) will always fit the data better 
than an Z max = 5 map (36 parameters) in the sense of having a larger value of p(D\a, M). 



p{M 1 \D) _ p(D\M 1 ) p(Mi) 
p{M 2 \D) ~ p{D\M 2 )p{M 2 ) 



(4.18) 




(4.19) 
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One can imagine, however, a stochastic background characterized by Vi m up to only / max = 5. 
In that case, an / max = 30 fit would introduce a great many unnecessary parameters, which 
means a much smaller value of 8a/ Aa offsetting the larger value of p(D\a, M). In addition, 
since the multipole moment estimators for different I and m are correlated with one another, 
choosing a large value of / max would have the undesirable effect of worsening the uncertainty 
associated with the Vi m up to Z max = 5. 

In the context of our search for an anisotropic stochastic gravitational-wave background, 
the model M is that a signal is present having multipole moments up to / max . The data are 
the measured cross-spectra Cft, and the likelihood function is given by Eq. (3.7), where the 
signal model is defined by multipole moments up to / max . Thus, the quantity we need to 
calculate is the marginalized likelihood 

p{{Cft}\lzaax) = 

J d{V lm }p({C ft }\{V lrn }, / max )p({P/ m }|/ max ) , (4.21) 

where p({Vim}\lmax) are the prior probability distributions for the multipole moments. Since 
the parameter space is large (a total of (/ max + l) 2 parameters), sophisticated Markov Chain 
Monte Carlo techniques (see [45]) may be required to numerically evaluate this integral. If 
the data turn out to be informative, then one has the much simpler expression 

PttC/dlUO « p({c f t}\{r lm }, u) n det A^" 1} ( 4 - 22 ) 

where Vi m are the maximum-likelihood estimators given by Eq. (3.10), and AVi m are char- 
acteristic widths of the prior distributions. Whether or not one can use this approximation 
depends on the actual data and the choice of priors. In practice, it may be possible to 
use limits from previous, less sensitive analyses to set the widths of the priors for (at least 
some of) the Vi m . In the absence of strong a priori knowledge, the widths of the priors will 
necessarily be large, reflecting our uncertainty in the values of the signal parameters. 

V. RESULTS OF SIMULATIONS 

In this section, we present the results produced by our data analysis code for simulated 
stochastic signals injected into simulated detector noise. We focus attention on analyses 
done in the spherical harmonics basis, as similar studies for the pixel-based decomposition 
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have already beeen discussed in detail in the context of the radiometer analysis [21, 32]. 
We find that the spherical harmonic analysis method can successfully recover simulated 
signals injected into simulated noise for several different types of stochastic gravitational- 
wave backgrounds, e.g., isotropic sources, dipole sources, point sources, diffuse sources, etc.. 
We also verify that the results of the standard isotropic and radiometer analyses are recovered 
as special limiting cases of the spherical harmonic decomposition analysis for Z max = and 
l max — »■ oo, respectively. 

A. Simulation details 

The simulations described in this section are made up of twenty-four jobs, each consisting 
of approximately one hour of data. Since the beam pattern matrix of the detector varies 
with local sidereal time, we chose the start time of each job so that the data are distributed 
(nearly) uniformly over a sidereal day. The twenty-four jobs are further broken down into 
one-minute segments (so r = 60 sec), on which the analysis described in Section III is then 
applied. 

The simulated time-series data are sums of simulated detector noise and simulated 
stochastic signals for several different angular distributions. The simulated detector noise 
are constructed so as to reproduce (on average) the design power spectral densities of the 
different detectors — in our case, the 4 km Hanford and Livingston LIGO interferometers (HI 
and LI) and the 3-km Virgo interferometer (VI). See Figure 5. The simulated stochastic sig- 
nals we consider include: no injection (i.e., just detector noise), an isotropic (i.e., monopole) 
source, a dipole source, a point source, two point sources, a diffuse source clustered around 
the galactic plane, and a diffuse source clustered around dec = 0°. Note that the dipole 
source is injected on top of a monopole of twice its amplitude, so that the signal power is 
positive everywhere on the sky. 

The spectral shape of the stochastic signal is taken to be constant (H(f) = 1) for all the 
injections and for all the analyses. The overall amplitude of the signals are different for the 
different injections, chosen to be large enough to be easily detectable in one sidereal day 
of total integration time. [52] Table I lists the expected values of Vqo/VAtt and ^ ? (^)| max ) 
either of which fix the scale of the various injections. The factor of 1/a/47t multiplying P o is 
included to allow direct comparison with the sky map plots of V(Cl) = X^Zm^Zm(^) shown 
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FIG. 5: The design power spectral densities used to simulate detector noise for the LIGO 4km 
interferometers (HI and LI) and the 3 km Virgo interferometer (VI) [46, 47]. 



Injection type 




V(tl)\ 

v ' 1 max 




(strain 2 /Hz/rad 2 ) 


(strain 2 /Hz/rad 2 ) 


Monopole 


5.6 x 10" 45 


5.6 x 10" 45 


Dipole 


1.1 x 1(T 44 


2.1 x 10" 44 


1 point source 


1.6 x 10~ 47 


4.1 x 10" 45 


2 point sources 


3.2 x 10" 47 


4.0 x 1(T 45 


Diffuse source (galactic) 


3.8 x 1(T 45 


2.0 x 10~ 44 


Diffuse source (dec = 0°) 


4.2 x 1(T 45 


2.0 x 10" 44 



TABLE I: Expected values of 'Poo/v / 4vr and 7 ? (^)| max f° r the different injections. 

later in this section, noting that Y 00 (£l) = 1/v^Lr. The maximum power values are given for 
easy comparison for the point source injections. 

The analysis code was then run on the simulated data, decomposing the relevant quan- 
tities with respect to the spherical harmonic basis as described in Section III. The main 
output of the analysis for a particular simulation are the spherical harmonic components of 
the dirty map X\ m and the beam matrix Ti m ^ m /. The maximum-likelihood estimates of the 
true multipole moments Vi m of the gravitational-wave sky are then obtained by inverting 
the beam matrix Ti rn ^ m i (either with or without SVD as discussed in Sec. IV B), and then 
applying that inverse to the X[ m to get the components of the clean map Vi m . (In what 
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Method 


Poo/v^ 


o"oo/\ / 4vr 




(strain 2 /Hz/ rad 2 ) 


(strain 2 /Hz/r ad 2 ) 


isotropic 


4.207339 x 10" 49 


3.209030411 x 10~ 48 


^max — 


4.207328 x 10~ 49 


3.209030408 x 10" 48 



TABLE II: A comparion of the maximum-likelihood estimates and error bars for the spherical 
harmonic decomposition code (Z max = 0) and the standard isotropic search. 

follows, a 'clean map' will mean the map constructed from the regularized inverse, unless 
we explicitly indicate otherwise.) 

B. Comparison with previous searches 

As an initial check of the analysis pipeline, we verified that the spherical harmonic de- 
composition code reproduced the results of the standard isotropic [23, 30, 31] and radiome- 
ter [21, 32] analyses in the limits Z max = and / max — > oo, respectively. For the isotropic 
comparison, we analyzed simulated data with no injection (just detector noise) and com- 
pared the Z max = results (i.e., the maximum-likelihood estimate Poo of the monopole 
moment, and the associated 1-sigma error bar ooo) to an identical analysis performed with 
the isotropic search code. The results, presented in Table II, show that the two methods 
give the same answers (to round-off error) for the isotropic component of the background. 

For the radiometer comparison, we analysed the same simulated data with no injection 
(just detector noise) with both the spherical harmonic decomposition code for different 
values of Z max , and compared the resultant dirty sky maps constructed from the Xi rn with 
the pixel-based map produced by the radiometer search code. Figure 6 shows that the 
spherical harmonic algorithm successfully reproduces the radiometer analysis in the limit 
of large Z max . For a radiometer pixelisation appropriate for the diffraction limited beam 
pattern at / ~ 1 kHz, Z max = 30 yields a good approximation. The difference between the 
^ma X = 30 map and radiometer map has fluctuations consistent with the angular scale set 
by Z max = 30 — i.e., the two analyses agree for angular scales accessible up to Z max = 30; they 
differ only for finer angular resolutions. 

In Figure 7 we show the sky map for this no-injection simulation that has been cleaned 
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FIG. 6: Left column from top to bottom: dirty maps of no injection (just detector noise) produced 
with spherical harmonics decomposition code for Z max = 5, / max = 10, Z max = 15, Z max = 20. Right 
column: dirty maps with ? max = 30, with the radiometer search code, the difference between the 
l miB = 30 map and the radiometer map. By reading top to bottom, one can see how the spherical 
harmonic dirty map approaches the radiometer map as /max increases. Residual fluctuations on 
the difference map appear consistent with the angular scale set by Z max = 30. 
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by the SVD algorithm, a SNR map for this clean map, and a histogram of the SNR map. 
It is readily apparent that the fluctuations are consistent with detector noise (no signal.) 

Thus, the spherical harmonic algorithm reproduces the isotropic and radiometer analyses 
as special cases, while allowing us the flexibility to consider more general cases, all within a 
single framework. 

C. Sky maps - Single baseline 

In this subsection we focus on simulations utilizing the Hl-Ll baseline. The sky maps 
constructed from the Vi m for the various stochastic simulations are shown in the following 
figures: 

1) In the top panel of Figure 8 we plot a dirty sky map for a point source injection 
with Z max = 20. While the location of the point source at (ra,dec) = (6hr, +45°) is readily 
apparent, the source is smeared and it is surrounded by artifacts arrising from the beam 
pattern function. If we attempt to produce a clean map by naively inverting T a p as in the 
second panel of Figure 8, we find that the 'clean' map is actually worse (less representative 
of the injection) than the dirty map due to singularities in the inverted matrix (as described 
in Section IV). In the third panel of Figure 8, we present a clean map derived using the SVD 
algorithm. The location of the point source is readily apparent, and the SVD algorithm has 
removed some of the artifacts and smearing associated with the dirty map. In the fourth 
panel we plot the associated SNR map. 

When comparing these maps one should bear in mind that the clean and dirty maps have 
different interpretations, and so the color scales have very different numerical ranges. In the 
illustrative case of Z max = 0, for example, X Q0 = Voo/ where X 00 is the dirty map and 
"Poo is the clean map. 

2) In Figure 9 we plot a clean sky map for an isotropic injection with Z max = 20. For 
this injection, V(Cl) = 5.6 x 10 -45 strain 2 /Hz/rad 2 , a value which is indicated by green in 
Figure 9. 

3) In Figure 10 we plot a clean sky map for a dipole injection with Z max = 20. 

4) In Figure 11 we plot a clean sky map for an injection of two point sources with 
/ = 20 

'max 

5) In the first panel of Figure 12 we plot an injection of a diffuse source clustered in the 
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FIG. 7: From top to bottom: A clean sky map for the case of no injection; next, the associated map 
of uncertainty, cr-p^y next, the associated SNR map; bottom, a histogram of SNR. The banded 
structure on the uncertainty map can be understood in terms of the directional sensitivity of two 
cross correlated interferometers. The rotation of the Earth ensures that the uncertainty is uniform 
in right ascension for a fixed declination. The blue histogram bars are data, the dark red line is 
a Gaussian fit (sigma=l, mean=l), the light green line is a maximum-likelihood fit (sigma=1.07, 
mean=0.06), and the dashed line is the 1 -sigma, error for 400 independent points. l m ax — 

20. 
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FIG. 8: From top to bottom: a dirty sky map of a point source injection at (ra, dec) = (6hr, +45°), 
a clean sky map (without SVD) of the same injection, a clean sky map (with SVD) of the same 
injection, a map of SNR. The source has maximum SNR of 49. Z max = 20. 

galactic plane. In the second panel we plot the clean sky map recovered from this injection 
using Z max = 20. 

6) One way to test that an injection is recovered successfully and without bias is to 
plot the injected signal map minus the recovered clean map. To do this, we must take into 
account that the clean map was produced using SVD (see Section IV B). That is, we need to 
compare the regularized extracted clean map V' a with a 'regularized' version of the injected 
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FIG. 9: Above is a clean sky map of an isotropic injection. Below is a histogram of SNR. The 
average SNR across the map is 9.1. Z max = 20. 
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FIG. 10: A clean map of a dipole injection oriented along the z axis. Z max = 20. 



map: 



Pi = r%^ 7 ? r (5.i) 

Here Y @ 1 is the Fisher matrix, (r /_1 ) a ^ is its regularized inverse, and V 1 is the injected map. 
In Figure 13 we plot V' a — V' a for the galactic injection depicted in Figure 12. 

7) In the first panel of Figure 14 we plot an injection of a diffuse source clustered about 
dec = 0° generated using Planck simulator [49] and HEALPix [48]. In the mid-left panel 
we plot the regularized version of this injection, and in the top-right panel, a clean sky map 
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FIG. 11: A clean sky map of two point sources, one at (ra, dec) = (6hr, +45°) (SNR = 81) and 
the other at (ra, dec) = (12 hr, -30°) (SNR = 76). Z max = 20. 
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FIG. 12: Above: a toy model injection corresponding to a map measured by the WMAP satel- 
lite [42] meant to mimic a diffuse source clustered in the galactic plane (6 = 0°). The map utilizes 
HEALPix [48] and the injection was simulated using the Planck Simulator [49]. Below: a clean 
map recovered from this injection. i max = 20. 

using Z max = 20. In the mid-right panel we plot V a — V' a . The apparent quadrupole moment 
visible in the mid-left and top-right panels illustrates the relatively low sensitity to I — 2 
moments using the Hl-Ll baseline. 
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sun 

FIG. 13: Top-left: the original injection. Mid-left: the regularized version of the injection V' a . 
Top-right: the recovered clean map P' a . Mid-right: V' a — V' a normalized by cr-pr^s- Bottom-left: a 
histogram of these residuals. The fluctuations in the residuals appear to be consistent with detector 
noise. Z max = 20. 

D. Sky maps - Multiple baselines 

Figure 15 shows the clean sky maps for a diffuse source distributed along the galactic 
plane (b = 0°) and Figure 14 shows the clean sky maps for a diffuse source distributed 
along dec = analyzed with single baselines (Hl-Ll), (HI- VI), and (LI- VI), and with the 
multi-baseline analysis (H1-L1-V1). 
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FIG. 14: Top-left: a toy model injection corresponding to a map measured by the WMAP satellite 
meant to mimic a diffuse source clustered around (dec = 0°). The map utilizes HEALPix [48] and 
the injection was simulated using the Planck Simulator [ L9]. Mid-left: a regularized version of the 
injection. Top right: a clean map recovered from this injection. Mid-right: the residuals V' a — V' a 
normalized by v-ptfiy Bottom-left: a histogram of these residuals. Note the apparent quadrupole 
moment present in the mid-left and top-right panels. This demonstrates the relatively low sensitity 
to I = 2 moments using the Hl-Ll baseline. l max = 20. 

VI. SUMMARY 

We have presented here a maximum-likelihood analysis method for estimating the angular 
distribution of power in an anisotropic stochastic gravitational- wave background. The basic 
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FIG. 15: Results from a multiple baseline simulation corresponding to the injection in the top 
panel of Figure 12. In the top panel is a clean map from the Hl-Ll baseline. Second from the top 
is a clean map from the Hl-Vl baseline. Third is a clean map from the Ll-Vl baseline. Fourth 
is a clean map produced from combining all three baselines (H1-L1-V1.) The final panel is the 
injected source. For all maps £ max = 20. 
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idea was to cross-correlate data from a network of two or more gravitational-wave detectors, 
exploiting time-of-arrival differences and the diurnal modulation due to the Earth's rotation. 
We derived maximum-likelihood estimators for the angular distribution of gravitational- wave 
power V(Cl) = ^2 a V a e a (Cl), decomposed with respect to any set of basis functions on the 
sky. We derived an expression for the beam pattern matrix T a p and discussed its relationship 
to the covariance matrix of the maximum-likelihood estimators V a . We described how sin- 
gular value decomposition can be used to regularize the inverse of T a/ 3, which was needed to 
remove the smearing effects of the beam pattern matrix on the measured ('dirty') sky maps 
X a . We also explained how the single-baseline (two-detector) cross-correlation analysis can 
be extended to a network of three or more detectors, thereby increasing our sensitivity to 
detecting a signal. In this paper, we focused attention on a decomposition with respect to 
a basis of spherical harmonics Y/ m (f2), for which the maximum-likelihood estimators V\ m 
represent the multipole moments of the gravitational-wave sky, and for which the standard 
isotropic and radiometer searches are recovered as special limiting cases. Finally, we illus- 
trated all these general results by analysing simulated data containing injected stochastic 
gravitational-wave backgrounds having different angular power distributions. 
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APPENDIX A: SPHERICAL HARMONICS 

Our convention for the spherical harmonics Yi m (9, <f>) follow [50]. Explicitly, 

= ^t 1 (f + m)! ir(c0s6,)e ^' (A1) 
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where P™(cos#) are the associated Legendre functions defined by 

( -\\m Al+m 

™ = y^(i-x 2 r /2 ^(x 2 -iy 



2 l l\ 



dx l 



The normalisation constants have been chosen so that 



and 



Note that 



and 



2/ + 1 {l-m)\ 



2ty pit 

d(j) / sin0d0 Yf m ,(e, (f>)Y lm (6, 0) = 6 vl 6 n 
o Jo 



Y 



l-m\P, 



:-irc(^ 



(A2) 
(A3) 



Y lm {-Cl) = Y lm (n -0,0 + vr) 
= (-l)^ m (0,0) 

= (-i) l Y lm (n). 

Expressions for the first few spherical harmonics (up to / = 2) are given below: 



Yao(e, 
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(A9) 




(A10) 



Y 



Y w {0A) 



-i\y, 




sm i 



7T 



3 

— cosy 

47T 




8tt 



sm i 



-up 



Y 



1 /15 



221 




4 V 2vr 



sm 



Y 



2U 



15 

— sin 9 cos Qe %i 
8n 



(All) 
(A12) 
(A13) 
(A14) 
(A15) 
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Y 20 (6, ( f ) ) = J^-( 3 -cos 2 e- 1 -) (A16) 



/ 1 Pi 

Y 2 -i(9,<l>) = \ — sin 9 cos 0e _i * (A17) 

V o7T 



F 2 ,„ 2 (^0) = y^sin 2 ^- 2 ^. (A18) 



APPENDIX B: USEFUL IDENTITIES 

The transformation property of the spherical harmonics (A6) and (A7) imply the following 
transformation property for the 7; m : 

iL(M = (-iy +m n,- m (f,t) (Bl) 

and 

7i m (-f,t) = (-l)\ m (f,t) (B2) 

= (-lTt-m(f,t). (B3) 

Similarly, the requirement that V(Cl) is real implies 

V; m = {-\TV^ m . (B4) 
Finally, Eqs. (Bl), (B2) and (B3), together with the definition (3.11) and (3.12) imply 

X: m = (-l) m AV m (B5) 

and 

Tlm,l>mi = for odd (/ + I') (B6) 

and 

(-i) m+m Fi - m j> _ m > = r* mtl , m , = r i / m / iJm , (B7) 

so j/ m / is Hermitian. 
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APPENDIX C: DETECTION STATISTIC 



In addition to estimating the individual components of an anisotropic stochastic back- 
ground, it is also possible to construct a statistic that is optimal for detecting the presence of 
a background characterized by a particular set of (normalized) angular components V a and 
spectral shape H(f). To show this, we note that in the presence of a signal, the components 
X a of the dirty map can be written in the form [32] 

x a = r af3 Vp + N a , (ci) 

where T a p and Vp are as before, and N a is an additive noise term composed of cross- 
correlated detector noise and stochastic signal components. In the weak-signal approxima- 
tion, the variance of the noise- noise cross-term, n\(f ,t)n 2 {f ,t), is much greater than that of 
the signal-noise cross terms, hl(f,t)n 2 (f,t) and h^f, t)ni(f, t), so to a good approximation 

- £ £ f.u, ') p l(/ ^ 2 ) (/ , t) f ««/. ')*></. o • «=> 

^ / 

Furthermore, when the detector noise is Gaussian and uncorrelated — an assumption that is 
well-approximated in practice — the N a are Gaussian-distributed with covariance matrix 

(N a N;)-(N a )(N;)^T Q p. (C3) 

To construct the detection statistic, we assume that the stochastic background has spec- 
tral shape H(f) and normalized angular components V a satisfying 

r a pV*Pp = i . (C4) 

The overall amplitude e of the background is given by V a = eP a . Then the probability 
density function for the X a in the presence of such a background is given by the likelihood 

1 



p({X a }\e) oc exp 



(Xp - eTpsP, 



8 



(C5) 



By the Neyman-Pearson criterion, the optimal detection statistic A is simply the maximum- 
likelihood estimator of e [51] — that is, 
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. (C6) 



€ = € 



The result, after a straightforward calculation is 



A = X a V* a , (C7) 

which has the form of a standard matched-filter. Note that the detection statistic A has 
zero mean and unit variance in the absence of a signal. In the presence of a signal whose 
parameters exactly match those of the signal model V a and H(f), the expectation value of 
the statistic is 

(A) = e. (C8) 

(The variance of the statistic is still unity in the weak-signal approximation.) In the special 
case of an isotropic background, A = Vqq/(Tqq, which is the signal-to-noise ratio for the 
standard isotropic search. 
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